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Abstract 

We consider a viscous incompressible fluid flow driven between two parallel plates by 
a constant pressure gradient. The flow is at a finite Reynolds number, with an 0(1) dis- 
turbance in the form of a travelling wave. A phase equation approach is used to discuss 
the evolution of slowly varying fully nonlinear two dimensional wavetrains. We consider 
uniform wavetrains in detail, showing that the development of a wavenumber perturba- 
tion is governed by Burgers equation in most cases. The wavenumber perturbation theory, 
constructed using the phase equation approach for a uniform wavetrain, is shown to be 
distinct from an amplitude perturbation expansion about the periodic flow. In fact we 
show that the amplitude equation contains only linear terms and is simply the heat equa- 
tion. We review, briefly, the well known dynamics of Burgers equation, which imply that 
both shock structures and finite time singularities of the wavenumber perturbation can 
occur with respect to the slow scales. Numerical computations have been performed to 
identify areas of the {wavenumber, Reynolds number, energy) neutral surface for which 
each of these possibilities can occur. We note that the evolution equations will break- 
down under certain circumstances, in particular for a weakly nonlinear secondary flow. 
Final!} we extend the theory to three dimensions and discuss the limit of a weak spanwise 
dependence for uniform wavetrains, showing that two functions are required to describe 
the evolution. These unknowns are a phase and a pressure function which satisfy a pair 
of linearly coupled partial differential equations. The results obtained from applying the 
same analysis to the fully three dimensional problem are included as an appendix. 
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1 Introduction 


Fully nonlinear travelling wave solutions in plane Poiseuille flow (PPF) have been discussed by 
many authors in recent years (Zahn et al [17], Herbert [9], Pugh k Saffman [14]) with emphasis 
p aced upon the secondary instability problem. This instability of the finite amplitude two- 
dimensional periodic flow to three-dimensional infinitesimal disturbances has been put forward 
traditionally as an explanation of transition. Orzsag k Patera [13] have also suggested that 
even at Reynolds numbers below the subcritical minimum for these nonlinear two dimensional 
solutions, the time scale for ultimate decay of the disturbance is sufficiently long for a small 
three dimensional perturbation to grow strongly. 

Subsequent work on superharmonic instability by Pugh k Saffman [14] has shown that a 
more complex structure is to be expected than a simple stability transition at the Reynolds 
number limit point of the neutral surface. They show that parameterization of the problem is 
important and that bifurcations to quasi-periodic flows exist at points on the upper branch of 
the neutral surface, leaving open the possibility of such finite amplitude solutions existing at a 
Reynolds number lower than the subcritical minimum for the periodic flow. 

It should be noted that other approaches have been discussed that do not rely on the 
existence of finite amplitude two-dimensional travelling waves. Although the possibility of 
resonant growth caused by linear mechanisms has been recognized for some time, it is only 
through recent work (Gustavsson [6] [5], Gustavsson k Hultgren [7], Butler k Farrell [2], Tre- 
fethen et al [15] ) concerning the initial value approach that the large amplification involved 
in three-dimensional problems has been revealed. The term “bypass” has been attached to 
these methods since they do not follow the more traditional idea of transition arising from 
two-dimensional Tollmien Schlichting waves with three dimensional effects appearing only at a 
secondary instability stage. 

In this discussion we return to the finite amplitude travelling wave solutions described at 
the beginning of this section, and subsequently develop an evolution equation for a phase 
instability of the (9(1) flow. Since these nonlinear solutions are used frequently in some areas of 
both theory and computation, it seems sensible to try and discover something of their stability 
and evolution. We show, in fact, that uniform wavetrains will not be observed under certain 
classes of initial condition as they are susceptible to slow scale effects, with the wavenumber 
developing both singularities and shock structures after a finite time in the slow scale. 

The method we use to determine the evolution equation is based upon the phase equation 
technique applied by numerous authors to Benard convection problems with (9(1) amplitudes 
(e.g. Newell , Passot k Lega [12]). These same methods have been applied to wave problems by 
Howard k Kopell [10], Whitham [16], subsequently applied explicitly to the Ginzburg- Landau 
equation by Bernoff [1] and utilized in an investigation of boundary layer instability theorv bv 
Hall [8], 3 ' 

A detailed discussion of the phase equation method is given by Hall [8], Bernoff [1] and is 
not repeated here. The essential idea is that we have a finite amplitude wavetrain solution, 
which is locally periodic in space and time, allowing wavenumber and frequency to be functions 
varying on appropriate slow scales. The resulting equations of motion can then be rewritten in 
terms of the new scales and a phase function that is related to the wavenumber and frequency, 
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through which it also satisfies a conservation equation. An expansion in terms of the slow scale 
parameter will now yield a leading order system that is a nonlinear eigenvalue problem. This 
relationship determines the local frequency of the wavetrain as a function of local wavenum- 
ber, Reynolds number (and indirectly amplitude) yielding both supercritical and subcntical 
equilibria for PPF. The next order problem will then provide a linear inhomogeneous system 
that determines the frequency correction term through a solvability condition; this technique 
also allows for continuation to higher orders. Now since the wavetrain evolves according to the 
phase conservation equation, we can, by expanding appropriately, give a slow scale asymptotic 
approximation to the evolution equation. 

In §2 of this discussion the above technique is applied explicitly to the finite Reynolds 
number two-dimensional PPF problem. In §2.1 we discuss the implications of the phase equation 
theory when applied to the stability of a uniform two-dimensional wavetrain. We show that 
the stability of a small wavenumber perturbation is governed by Burgers equation, 

A t + AA C = ± A^ , (!) 

for 0(1) problems that correspond to distinct points on either the upper or lower branch of the 
neutral surface and away from the linear neutral curve. In §3 we discuss the stability of the 0(1 ) 
flow to an amplitude perturbation, showing that, for these length scales, nonlinear terms are 
not introduced and that the amplitude equation is simply the heat transfer equation. Section 4 
provides a short description of the numerical methods involved in the uniform stability calcula- 
tions, presenting the 0(1) results which are consistent with those of other authors Herbert [9], 
Pugh &: Saffman [14]. Computational results, which show the behaviour of the viscous diffu- 
sion term from Burgers equation, are presented for differing leading order problems. Section 5 
returns to the analysis of a uniform wavetrain and briefly discusses how the theory breaks down 
for weakly nonlinear secondary flow solutions. In §6 we redevelop the phase equation theory 
for a three-dimensional problem and consider the stability of wavenumber perturbations m 
the limit of a weak spanwise dependence, finally in §7 we discuss the implications and future 

extensions of the work. 


2 Formulation of the phase equation approach 


We wish to consider a finite amplitude solution to the plane Poiseuille flow (PPF) problem, in 
the form of a travelling wave, then allow for a slow modulation on the new scales, 

X = 8x and T = St . (2) 


We now use the methods presented in Howard k Kopell [10], Bernoff [1] and applied to asymp- 
totic suction boundary layer flow by Hall [8]. This analysis follows closely the finite Reynolds 
number case of Hall [8], except for a few technical differences associated with bounded domains, 
namely the need for a pressure eigenfunction term (discussed later). We first introduce a phase 
function, 0(X, T) = 86{x,t), which allows a definition of the local frequency and wavenumber 


as 




(3) 
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where a - a(X,T) and ft - f t(X,T) are allowed to vary on the slow scales. Thus the partial 
derivatives transform as 


— _d r d d_ 

dx °de + dt 



c d 

+ S dT ’ 


and a conservation of phase condition must be satisfied, 


(4) 


da dfl 

dT + dX = °* ( 5 ) 

We can now develop a perturbation scheme about the fully nonlinear leading order solution 
by introducing a slow scale expansion of the stream function, 


0 = 0o + Sipi + . . . , 

which forces a similar expansion for the frequency 


(6) 


ft — flo T <5X2} 5 


(7) 


where 0o - 0 + 0o the basic flow potential plus a leading order spatially periodic flow. We 
shall use the vorticity equation formulation to describe the flow, 


g(V 2 0) d(V 2 0, 0) 

dt d(x,y ) 


i v4 0 = 0 ’ 


( 8 ) 


with a Reynolds number, Re, defined as where U* is the centerline velocity, 2 h the plate 
separation and !/ the kinematic viscosity. Thus at leading order, using the above expansions 

ohf«si« ' 


-S7 O V 2 0 Off + a0 O!/ V 2 0 Off - Q;0 O $ V 2 0 Oy - -^-V 4 0o = 0 , 


Re 


(9) 


with subscripts {0, y,...j denoting the respective derivatives, where it is unambiguous to do 
so, and 


d 2 d 2 


v 2 = a 2 ^- + 


(10) 


d9 2 ' dy 2 ' 

The boundary conditions are simply impermeability and no-slip at the parallel boundaries. In 
terms of the streamfunction (at this order) this gives 


000 = 0Oy = 0 at y = ±1 . (ii ) 

Similar situations have been discussed previously by Herbert [9], Pugh & Saffman [14], and this 

nonlinear eigenvalue problem can be solved using the same numerical techniques, here we have 
an eigenrelation 

= fi 0 (aq Re) , (12) 

which defines a “neutral surface” in {cv, Re, Amplitude} parameter space. 
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Now at next order, after some rearrangement, we obtain 

r „ * dSl 0 dV 2 i/> 0 djfo, 2aipoe a + rj>oe) 

L{rpi}= Slo&cxxpoee* + $w) + q q g^r~ + <* d{0,y) 

+ + J_(4aV 2 ^a + 4a 2 ^oeee + 2V 2 tM] V 2 ^oe 


( 13 ) 


d(a,y) ' ^ 

with boundary conditions 0 ly = r/>i* = 0 at y = ±1. To obtain the above form, (13), we have 
introduced the operator L defined as 

d ^ c) \ 

L = -^oV 2 ^ + axP 0y V 2 ^ + <*V 2 i/>o^ - "^V 2 ^ - «V 2 t/» 0y ^ - -^V 4 
and the d/dT term has been replaced using the conservation of phase as 

d 


(14) 


d . da d 

dT '-^(^0 + ^,+...) 


dX da 


(15) 


We now consider the mean flow, which from the Navier-Stokes is governed by 


1 d 2 


Uoo 


( d u o j , (c)duo£\ 

= E(-'»^ + %^-j + 


Re dy 2 


du o 
dy 


dq-i 

dX 


Here we have expanded the velocity field and pressure as 

u = Ug(X, T, 9, y) + <$Ui (X,T, 6,y) + . . • , 
p = [GX + T)]«-' + \po{X,T,0,y) + qo(X,T)} + 


with 


Uo = (0, of + Uoo + £{noi eii# + VLo1 e *'*} ’ 2oi = ( U( 0’ v °^ 

j=l 


(16) 


(17) 

(18) 

(19) 


where U denotes the non-dimensionalized basic flow, G - the basic driving pressure 

gradient and (c) a complex conjugation. At next order we obtain 


1 d 2 u w dpoo , dQoduoo pdu 0Q 

Te -W = ax + aT -^ ax + ' u ' w + dy v '° 


oo , 

+ E ( “") 


dutf duV . SufcK dqo 


J = -CO 


dX 


0 j , vu, lj ~“0j 

+ v 0j a.. T Ul i' 


dy 


dy 


+ 




( 20 ) 


where poo is determined from the 0(1) problem and 


u, = u,„ + + uti’e-' 9 } , U>. = (««. • 


( 21 ) 
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the summation term in (20) is performed with negative subscripts denoting a complex conju- 
gation. The equation of continuity at this order, 0(5), is 


du x du 0 dvx . 

“¥ + ax + l^~ 0 ' 

which, when considering the mean flow terms only, reduces to 


dvxo duoo 

~d^ = ~~dX' 


( 22 ) 


(23) 


a first order equation for the mean flow correction to v, required to satisfy the impermeability 
conditions at both walls. This difficulty was anticipated earlier and is resolved by the introduc- 
tion of a further, slow scale dependent, pressure expansion 


q- 1 (X,T)8~ 1 +q 0 (X,T) + ..., (24) 

producing the extra term dq-x/dX in (16) which is chosen to satisfy (23). Thus, for imperme- 
ability at both boundaries, we must satisfy 

r+ J duoo 
/-i ~dX 


i: 


dy = 0, 


(25) 


which fixes the streamwise flux through the channel and hence determines q_\ as a function of 
X at given T. This is equivalent to solving the vorticity equation for the mean flow correction 
with the same boundary conditions rp 0 o = dtpoo /dy = 0 at y = ±1, where tp 0 o is the part of the 
stream function having zero mean with respect to the phase, 8. 

Now we have an 0(5°) problem that can be solved numerically to give fio(c*, Re) and an 
expansion in terms of 6 giving a further system (13), which is used to compute flj at a given 
neutral surface point. The homogeneous form of (13) is solved by ipoe (corresponding to the 
existence of a translationally invariant solution, since any arbitrary constant may be added to 
the phase) and so fly is determined by a solvability condition at 0(5). The form of (13) is 

O r\ 

Q; ^( /1 3) + ^( B a) + (:7 3 = H’ ( 26 ) 

so premultiplying by the adjoint vector, r r = (rq, . . . ,r 6 ), and integrating by parts gives 

-o^(2l T r)-| ; (B T r) + C T r = 0 (27) 

as the adjoint equation for the homogeneous form of (26), with 

q = (^ 1 , <*019, 01 y, V 2 01 , <*V 2 V> 10 , V 2 0 !j,) T , (28) 

000 0 O' 

0 0 0 0 0 

1 0 0 0 0 

0 0 1 0 0 ’ 

0 0 0 0 0 

0 0 0 - 4-0 



1 

0 


0 

0 
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and 


C = 




■ 0 

0 

0 

0 

0 

0 ' 




1 

0 

0 

0 

0 

0 



B = 

0 

0 

1 

0 

0 

0 



0 

0 

0 

0 

0 

0 




0 

0 

0 

1 

0 

0 




. 0 

0 

0 

0 

0 

1 

Re 


0 

-1 


0 


0 


0 

0 

0 

0 


-1 


0 


0 

0 

0 

0 


0 


-1 

0 

0 

0 

0 


0 


0 


-1 

0 

0 

0 


0 


0 


0 

— ; 


0 -VVoy otW 2 ip06 0 (t^oy - -Oc^oe 


( 30 ) 


(31) 


The boundary conditions are also determined by the above process and are easily shown to be 

r 5 = r 6 = 0 at y = ±1 , (32) 

plus a periodicity condition. Thus the solvability condition is 

r 2 it y +1 


r 2 ir r + 1 

/ / r 6 H 6 dy dO = 0 , 

Je = 0 Jy=- 1 


(33) 


giving _ 

= f +1 M(0,y)dyde}{ T f +1 -V^wdyd*} , (34) 

OX l J6=0 Jy=-\ ) Jy=-\ > 

where M is the term in the square brackets in (13) and H = (0, . . . , 0, H$) T . 

Given Hi we can now write down the phase conservation condition upto 0(8) in the form 


da 

3T + 


da 

IfadX 




(35) 


This is an evolution equation for a = a(X, T), which we could in principle continue to any 
order through this expansion scheme; rather than solving this directly we shall restrict our 
attention to a somewhat simpler problem. 


2.1 Stability of Uniform Wavetrains 

Here we wish to discuss the stability of a fully nonlinear uniform wavetrain to a small, slow 
scale, wavenumber perturbation. Hence, following the above expansion scheme, we obtain an 
0(6°) problem with a phase function of leading order form 

0 O = aox - f2o(c*o)2 , (36) 

where {<* 0 , Re) defines a point upon the neutral surface at 0(1) disturbance energy, with 
associated frequency fl 0 - If we now perturb the wavenumber of this uniform solution with a 
slowly varying function A, 

a = c* 0 + A(X, T ) , (37) 
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and perform the above analysis, the evolution equation reduces to 


-d A / V dQo d A / \r rp\ , nO 9^2 

dT A[X,T)+- s -—A(X,T) + S w+ 6' w + ... = 0. 


( 38 ) 


Now we introduce a Taylor series expansion for SI, about the uniform solution, and a transfor- 
mation of the stream wise coordinate to a new frame of reference, moving with speed fif(ao) 
to obtain uv n 


dA 

V— + 


a 2 n 0 (oo),3A a 2 A / sa 

ar T da 2 ap = af. 


,.d ‘ A yji’i 

“ae •*-*) ■ 


(39) 


Here we have defined r. £ and <&(a) by 


so that 


T = T]T, 

t = X-n’ 0 (a 0 )T, 

n. = &#m. 


£ 


- 

d_ 

d(' 


(40) 

(41) 

(42) 


(43) 

(44) 


We therefore obtain a leading order balance for rj = 6 and A ~ 0(8). 

Now an appropriate rescaling will reduce this evolution equation to the standard form for 
Burgers equation, namely 


A t + AA { = fxA a , (45) 

which has well known properties. This approach again follows that given by both Hall [8] 
and Bernoff [1], where Burgers equation has been derived by applying this technique to both 
the asymptotic suction boundary layer and the Ginzburg-Landau equation. Whitham [16] has 

discussed the dynamics of this equation in detail, it has an exact solution via the Cole-Hopf 
transformation, 


A = -2;i— log(^) , (46) 

which removes the nonlinear term reducing the evolution equation to the form of the heat 
diffusion equation, 


</v = nil’# + C(t)i/j, (47) 

where C(t ) is set to be identically zero since this merely corresponds to a scaling of the depen- 
dent variable. Now for the initial value problem, with known A(£,0) = F(£), it is possible to 
obtain the analytical solution 


T ) = 


1 

y/XjFJvr 



(j - riY 

A (IT 




(48) 
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through an application of Laplace transform methods. Therefore the solution is 
A«,r) = {£ 


r+°° t-v c - D/2 ^ 

r 


where 


D = f F(r))dri + 
Jo 


K/r 

(49) 

U - v? 

(50) 


2t 


This can support both shock structures and singularities at finite times. For a positive 
diffusive term on the right of (45) we have a bounded solution which will decay for a local- 
ized/periodic disturbance. For a negative right hand side the solution is diffusively unstable and 
will become singular at finite time, indicating that the slow variation assumption is no longer 
appropriate and a return to the full equations of motion is required. Weak shock structures are 
discussed by Bernoff [1] who notes that for a small monotonic wavenumber variation such that 

(51) 


and 


lim a = ot- 

X— + — OO 


lim a = a+ , 

r— *■+00 


(52) 


(53) 


then Burgers equation applies as a leading order form for the evolution equation if 

Aa = a+ — a_ <C 1 * 

Here we also require that the unmodulated wavetrain corresponds to a distinct point on either 
the upper or lower branch of the neutral surface at 0(1) amplitude. Now this wavenumber 
variation will eventually become concentrated into a weak shock structure of width 0( A J in 
the unsealed streamwise coordinate, if 

(54) 


SV 0 (a 0 )(a+ - «-) < °, 


moving with speed 


c = 


fl(a + ) — fl(a_) 


Aa 


(55) 


a discretized form of the group velocity. If this variation in wavenumber (Aa) increases then 
the slow scale assumptions are eventually lost (as with the finite time singularity case) and the 
evolution of the wave system is governed by the full equations, namely the Navier Stokes. 

In this case (stability of uniform wavetrains) we should also note that the leading or er 
problem is simplified since there is no slow scale dependency for the wavenumber Thus the 
effect of the extra pressure term dq-i/dX (now a constant) is to induce an additional parabolic 
velocity profile into the mean flow correction and corresponds to a scaling of the Reynolds 
number at fixed amplitude/ wavenumber. So the condition (25), which fixes the flux through 
the channel, effectively determines a unique parameterization of the problem (as discussed by 
Pugh k Saffman [14]) in a self-consistent manner. 
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3 An Amplitude Perturbation Approach 

We now show how the phase equation method described previously is distinct from a more 
typical amplitude perturbation approach. In this method we solve the same leading order 
problem for the uniform wavetrain, but it is now perturbed by an eigenfunction with a slowly 
varying amplitude B = B(X , T). Again we introduce the slow scales 


* 

II 

o-> 

R 

(56) 

and 


T = Sty 

(57) 

together with 


II 

1 

(58) 

a new moving coordinate system, and a further timescale 



T = ST. (59) 

This slower timescale is known from the previous section but could otherwise be determined 
from the final solvability condition. Here c g is a group velocity and given values for a and ft we 
can expand in terms of 6 and a phase variable 0 — ax — fit. Now, seeking a solution analogous 
to that in Hall [8], we expand the stream function as 


where 


+ Stpi + . . . , 


V>1 = B(X,T)a 

with vp corresponding to the basic parallel flow, and 


dtpo 
30 ’ 


_ [B(X,T)afd^ o , 7 

V ” n! ’ n - 2 ' 


(60) 

(61) 

(62) 

(63) 


We now return to the vorticity equation (8) and substitute the above expansions along with 


d 

d d 


dx 

a d9 + ~dX' 

(64) 

d 

r.8 r 9 8 


dt 

— S' c g ~ * T S — r , 

do 9 dX dT 

(65) 

A? 

d 2 2 d 2 


V 2 = 

(X 

d y 2 ^ de* ’ 

(66) 


to give a leading order the form of the vorticity equation 

0(6 ) : — f2V 2 0o0 + — oc'tpoe^'^oy — — 0 . 

Re 


(67) 


9 



At next order 


0(5) : -nv^is + a0 Ov V 2 0i0 + onl>iyV 2 xj>M 

— a0o0V 2 0ij, ~ <*0i0V 2 0o y — V 4 0i = ® ’ (®8) 

with solution (62), which is an amplitude perturbation of the underlying periodic flow; in this 
formulation 0 O is independent of the slow scales X,T. The required group velocity is now 
determined from the next order system, which can be rearranged more clearly as 


Lo(6) {02} — Cc 


2ot£lipoeee + c 3 V 2 0 O 0 — 0oj,(2a 2 0 O 000 + V 2 ^ 0 fl) 


4a * 


+ xi)oe{2a 1 xpoee y + V 2 0 O3 ,) + — V 2 ^o^ 


B 


X ’ 


(69) 


with the L 0 (s) operator defined by the 0(6) equation and once terms proportional to B 2 have 
been eliminated by taking djdO of (67). We observe that 0 2 = aB^'tp oa and as expected c g 
corresponds to 50 0 (a)/5a in the phase equation approach, this follows by taking d/da of (67) . 
Note that an additional multiple of the homogeneous solution to 0 2 will not alter the solvability 
condition at next order, but will contribute to an amplitude equation at higher order. 

We also must remember that the additional pressure term (<7_i, discussed in §2) is still 
required in the leading order mean flow problem. Obviously a similar pressure term, dq-i/da , 
is necessary at 0(6 2 ) but it is not until 0(5 3 ) that the condition determining 1 is obtained. 
Now continuity of mass, at 0(5 3 ), requires the same constant mass flow condition to be satisfied 
for impermeability of the boundaries. Thus we determine q- 1 , appearing in the leading order 
problem, in the same manner as discussed in the phase equation approach. 

The same process can be repeated for the problem at next order which becomes 

L 0( sM 3 } = BBx{K(9,y)} + B f {-aV 2 t/> 0 *} + £**{m( 0, y)a} , (70) 


after again eliminating the B 3 terms by taking 5 2 /50 2 of (67). The M expression, in (70), is as 
defined previously in the phase equation analysis with {f} 0 ,5 a ft 0 } replaced by {fi,c s } and K 
is given by 


K = 2 a 3 niPome + c 3 a 2 V 2 0 o *0 + 2 a 4 

2 C>(0o, V 2 rp 0 e) , 3 5(000, V 2 0 Oa ) , 0 4 5(0 O 0,0 000 ) 

+a —a(¥^r +a 8(m +2a ^^r 

2 m o,,vvo) 

+ 9(0, ,) +a 9(0, y) Re 


Although at first sight (70) appears to be Burgers equation, reproduced through an ampli- 
tude perturbation approach, it does in fact reduce to a simpler form once we have observed that 
a particular solution is available to remove the nonlinear term from the solvability condition; 


'O(S) 


{a- 


50 Oo 

de 


} 


= K. 


(72) 
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Similarly by an inductive method we can show that a rescaling does not reintroduce nonlinear 
terms (which will be of the form B n dB/dX) into the amplitude equation, since we can develop 
a general particular solution; 

d n il> 0a 

n! dd n ' ( 73 ) 

If we now follow the same method outlined in section two we can obtain a solvability 
condition at this order, namely 

Bf J Jq — d# dy + Bjfjjr j J Mr§ d 0 d y = 0 . (74) 

This is the heat transfer equation^ the solution is characterized by the sign of the diffusive term, 
given by 

{£-J~ Mr,ieiv }{C-i £ < 75 ) 

in the previous notation. So we note that these amplitude perturbations are of less interest 
than the corresponding phase instabilities; since they are governed by the heat equation the 
solutions will simply decay exponentially to zero or grow becoming singular in a finite time. 
Thus we have stable decaying solutions for $(a) < 0 but when $(o) > 0 we must at some stage 
return to the full evolution equations to determine the development as higher order spatial and 
temporal derivatives are reintroduced. 

4 Numerical Methods 

4.1 The leading order problem 

In solving the leading order form of the vorticity equation we look for solutions that expand as 

°o 

V>0 = $ + Y V’one*"* , (76) 

n——oo 

where ip 0 is the leading order term in the expansion of the stream function ip, (6), and 0 = DV 
• The vorticity equation (8) reduces to, after substitution of the above form, to 

4-(Z) 2 - n 2 a 2 ) 2 ip 0n - tan | ^ {D 2 - n 2 cx 2 )xp 0n - ip 0n D 2 U J 

°o , 

+ ia Y 1 ( n ~ "i)V> 0 n-m(£> 3 “ rn 2 a 2 D)ip 0m 

m= — OO ^ 

-mDip 0n -m(D 2 -m 2 a 2 )0om) = 0, (77) 

for n = 0, ±1, ±2, . . . with boundary conditions ip 0n = Dtp 0 „ = 0; where D = d/dy. Now we 
solve this with a truncation of the Fourier modes and a Chebyshev expansion in the y-direction 
for each harmonic, 

Jpon(y) = — ^ + Y a nrT r {y), (78) 

* r=l 
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Figure 1: A cross section of the neutral surface at a = 1.1 N h - 1,2. 


along with U = \{ 1 - T 2 )■ This formulation has been applied previously by Herbert [9] to plane 
Poiseuille flow with a constant pressure gradient rather than the constant mass flux condition 
we apply in this case. 

If the harmonics are truncated at N h , and the Chebyshev series at N c , this yields a numerical 
problem that can be solved using Lanczo’s r-method with Nh.(N c + 7 ) nonlinear equations plus 
a coupled mean flow problem. We also note that the computational task can be simplified by 
assuming the symmetry condition 

Ip n(y ) = (-1 ) n+1 *Pn{-y) , ( 79 ) 

together with requiring a real solution, 

rp-n = ip+l . (80) 

The r-method is essentially equivalent to determining the higher Chebyshev coefficients through 
the boundary conditions, however we can only replace two of the dynamical equations with 
boundary conditions and therefore have to retain at least one r element. Various methods were 
investigated for the solution of these simultaneous nonlinear equations and finally a Newton 
iteration technique was found to provide the best convergence over large amplitude ranges. We 
write the nonlinear system of equations as 

f(x) = 0, (81) 
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Figure 2: Variation of frequency with E at a = 1.1 . 


where 


X — (Do? a ll? Ci2, . . . , 0*1 Nci a 20, ^21, - • • , 0'2N C ’> a N h 0, • • • , a N h N c ) 7 


and at each iteration level solve 

■fc = -£(&) , 


to give 


2&b+i = + z , 

where J is the Jacobian 


M 

dh 

Ih. 1 

dxi 

dh 

M ‘ 

fe 

dx\ 

dx 2 

dXm 

dim 

dim 

dfm 

_ dx\ 

dx2 

dx-m 


( 82 ) 

( 83 ) 

( 84 ) 


( 85 ) 


The number of unknowns in this solution procedure can be effectively halved by applying the 
symmetry assumptions discussed above. From this nonlinear eigenvalue problem we determine 
the relationship between fi 0 and {a, Re}, to do this we specify an amplitude of disturbance by 
choosing a value for the first Chebyshev coefficient, ai 0 , in the expansion of the fundamental 
mode. Without loss of generality we can assume that this amplitude measure is real since this 
corresponds to a unique determination of the phase. Here the initial guess for a solution of the 
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Figure 3: A cross section at Re = 3500 for N h = 1,2. 


system, at a general neutral surface point, is derived from an interpolation of previous “nearby” 
solutions. Since a value for one of the Chebyshev coefficients is specified we must replace it by 
some other unknown in the Newton iteration scheme, typically the frequency. 

As discussed in §2 we must also iterate upon the additional pressure constant, dq~ \fdx in 
the mean flow equation (16), to satisfy the constant mass flux condition. At each level of the 
iteration scheme described above we compute x fc , then the corresponding mean flow problem 
is solved directly, including the constant dq-i/dx, from (16). We continue the iteration, as 
described above, until some measure of convergence is satisfied, which in this case we choose 
to be that the frequency correction and largest change in Chebyshev coefficients are less than 
preset tolerances. 

In figure 1 we show a typical cross section of the neutral surface at a fixed wavenumber, for 
which there are no linearly unstable modes. These results are consistent with those obtained 
by Pugh & Saffman [14]. The disturbance energy here is defined as 

N h 

£=££n, (86) 

n=l 

where 

E n = Y I-i (l' D ^ 0n l 2 + l m *V’°n| 2 )dr/, (87) 

a normalized energy measure for each harmonic. 
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Figure 4: The neutral surface at a = 0.95 


In figure 2 we show the variation of the frequency for the same cross section, and figure 3 
displays the form of the cross section for a fixed subcritical Reynolds number. Finally in figures 
4 and 5 we show several cross sections for separate wavenumber choices of {0.95, 1.08} showine 
the effect of higher harmonics. 5 

We found that N c = 40 was generally sufficient to give converged results for a chosen number 
of harmonics in this Reynolds number range. Although the inclusion of higher harmonics still 
altered the mid-range amplitude results the qualitative behaviour remained the same, as noted 
in Herbert’s investigation [9] of the slightly different problem of a constant pressure gradient 


4.2 The 0 ( 8 ) problem 

To determine the sign of the diffusive coefficient on the right of Burgers equation we need to 
solve the next order problem. In §2.1 the phase equation analysis of the uniform wavetrain 
problem showed that the frequency correction term was determined by the solvability condition 


Q 


da 

dX 



"2 y 1 ^ f y*^ ?r r+i r „ 

-V tpoe r 6 dy d# > l / fi 0 (2a^w« + j>w) 

) d t/~0 Jy— — l L 
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<9V 2 V>0 d(t/>o,2c»/>oga + V>Qg) dijkh V 2 ^o) 

+ O 0o , ^ +a ^ y ) 9(a,y) 

+ Re~ l {4aV 2 j>o« a + 4a 2 V><*w + 2V 2 ^}]r 6 dt/dflj . 


Here V> 0 is the basic parallel flow plus two dimensional periodic flow and r 6 solves the adjoint 
problem, which can be shown to be 

— Re~ l V 4 rg — {orrf>oy ~ ^o)^ 2r 60 + r 6 y 

-2a^(tt 2 ^ - ]^) r s + 2ar6B y{ Q2 ~iw ~ d^)^° = °’ ^ 

with boundary conditions r 6 = r w = 0 at y = ±1. However in our calculations we do not solve 
the linear homogeneous problem (89) but rather define a new composite nonlinear problem 
that is an inhomogeneous form of the order one system, this can be computed directly after the 
leading order solution using the same routine. We solve 

N{j> o; Ho} = <5M(0, y) + 0(8 2 ) , (90) 

where the operator N and the new parameters are given by 

N = £o(i) + 8Lq(s) i 
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Reynolds Number, Re 

Figure 6: The variation of <f> at a = 0.95 &; = 2. 

fio = flo + <5 fix , (92) 

V’o = V’o + Sfa ) (93) 

with M and To(.) as defined previously. Hence after solving the 0(1) problem we compute 
the inhomogeneous terms produced by M then iterate on the amplitude of the new composite 
solution at fixed {a, Re }. The frequency correction term is then determined by an application 
of the condition fix = dCl 0 /d6 at 8 = 0. 

A point to note is that we divide the neutral surface into two sections which are referred 

to as the upper and lower branches. We can see from (88) that we should expect a singular 

behaviour for the coefficient of the diffusive term as we approach two other regions, namely 
where the lower branch connects with the E = 0 plane and where the upper and lower branches 
join at finite amplitude. The singularities are due to the way in which the amplitude and 
frequency vary with a at these points and further work is required near these regions. 

We present in figures 6 and 7 the behaviour of the diffusion coefficient in the evolution 
equation (as the 0(1) solution varies along the lower branches shown in figures 4 and 5) since 
this essentially characterizes the solution; we plot <J>(<* 0 ) where fix = 4>(c*)A*. Since the linear 
neutral curve exists for each of these wavenumbers we will find a singular behaviour of $ at 
a higher Reynolds number where the neutral surface intersects the E = 0 plane as discussed 
above. 

We have carried out calculations (of <b), w r ith more harmonics, for a number of typical 
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Figure 7: The variation of $ at a = 1.08. 


parameter values and find that such results have the same qualitative behaviour. In fact the 
main effect of the higher harmonics is to alter the cross section shape as shown in figure 3 
thus altering the position of the singularities when evaluating $ for varying Reynolds number, 
wavenumber or amplitude. In figures 8 and 9 we show the behaviour of the coefficient as a 
function of the wavenumber a, the vertical lines in figure 8 show the critical wavenumbers 
predicted by linear theory. Finally figures 10 and 11 show the behaviour for solutions that 
correspond to points upon the upper branch of the neutral surface at a. = 1.1, Nh = 1 und 
O' = 0.9707, Nh = 1,2 respectively. 


5 The Evolution Equation as E — > 0 

As noted in the previous section the phase equation method is only valid for finite amplitude 
periodic states; if we allow the leading order solution to approach the linear limit then the 
coefficient 4>, in Burgers equation, becomes singular. A simplistic approach would be to return 
to the uniform wavetrain problem and allow the 0(£°) system to be a small amplitude solution 
with a phase function of leading order form 

— &o x — fIo( a o)f - (94) 
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Figure 8: The variation of $ for constant supercritical Re with N h = 1 . 


In this formulation {c*o, Re} is a point on the lower branch of the neutral surface in the weakly 
nonlinear region with 

<*o = <*„ + «, |a|<l, (95) 

where a n is a wavenumber corresponding to a point on the linear neutral curve. If we perturb 
this solution with a slowly varying function A, 


a = a 0 + A(X,T), (96) 

then, as we approach the linear neutral curve, we know that the neutral surface has the scalings 

A = 0( |a-a n |s), (97) 

fi = fl n + 0(a — a n ) , (98) 

for some amplitude measure A and where 0 n is a frequency associated with the linear solution. 
The full form of the evolution equation is 


fr + (p°( a °) + Afi oM + •••)!§ = + ^ 2 ( a ) + •••}> (99) 

or equivalently, 


*7 q t + ^o( a o)A— + |fIo(a 0 )A 2 — + . . . j = £<I>(a;o) 


d 2 A 

ae 
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Figure 9: The variation of $ for constant subcritical Re. 

+ *{$'(<*<>) + $'(a 0 )A^| + ■••}+ + • • • ’ ( 10 °) 

after applying the substitutions (40)-(42). 

Therefore approaching the linear neutral curve may reintroduce other terms from the ex- 
pression 

ggi _ (94> 
d£ da 0 

since applying the scalings (97)— (98) to equation (88) we see that 

$(<*„) ( 102 ) 

a 

and the above terms can be of the same order of magnitude when 

jj ~ ~ (5 2 ~ Ct . (103) 

However the theory in this limit is more complex, as we approach the linear neutral curve the 
amplitude will no longer be determined explicitly from the leading order eigenvalue problem. 
So instead of assuming the weakly nonlinear scalings (97)-(98) we must determine the real 
amplitude of the weakly nonlinear uniform wavetrain by solving an appropriate nonlinear partial 


dA 


d$ . d 2 A 




dt 


da 0 di 2 




( 101 ) 
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Figure 10: The variation of $ for upper branch solutions at a = 1.1. 


differential equation. This amplitude equation will in general be coupled to a limiting form of 
the phase equation. Thus as we move from a fully nonlinear to a weakly nonlinear uniform 
wavetrain solution (for the leading order problem) we obtain two eigenfunctions for the problem 
and so must satisfy two, coupled, real partial differential equations (for the two-dimensional 
case). These two real equations may then be combined to recover a single complex evolution 
equation equivalent to that obtained by an application of the Stuart— Watson approach to the 
plane Poiseuille flow problem. Furthermore, in this limit, we no longer have a small parameter, 
6, that is arbitrary, it will now be related to the magnitude of the wavenumber displacement 
from the linear neutral value, a n . 

The connection of the phase equation method to the appropriate weakly nonlinear theory 
proves difficult in general, indeed the analysis still remains incomplete for convection problems in 
which phase-equation methods have been applied for some years. The recent works of Newell, 
Passot & Lega [12], Newell Sz Passot [11] and Cross Sz Newell [3] have discussed the points 
above and gone some way towards resolving the difficulties, although they restrict attention 
to much simplified model equations for convection (e.g the Swift-Hohenburg equation). In 
particular Cross and Newell [3] show how the evolution equation for such a model system can 
be matched with the Newell- Whitehead-Segel (NWS) equation in the weakly nonlinear limit. 
In their analysis they demonstrate that the limiting forms of the phase evolution equation and 
the amplitude equation (required in this limit) form the imaginary and real parts of the NWS 
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Reynolds Number, Re 

Figure 11: The variation of $ for upper branch solutions at ex = 0.9707. 
evolution equation respectively, hence they can be combined to match with previous results. 


6 A 3-D Phase Equation Theory 

IVe now apply a generalized (following Howard and Kopell [10]) form of the two dimensional 
phase equation method to a spanwise dependent problem. To achieve this we first introduce a 
further spanwise scale 

Z = Sz , (104) 

and generalize the phase as a function of three variables, 0(x, z , t) = 0(^, Z , T)/S, now defining 
a spanwise wavenumber, /?, by 

/?=©*. (105) 

The consistency conditions in this approach are 

0:7- + fix = 0 , (106) 

c*z — Px = 0 , (107) 

/3 r + ft z = 0, (108) 

where both wavenumbers and frequency are functions of the slow scales X , Z , and T . We restrict 

attention to the stability of uniform wavetrains, considering a leading order problem that is an 
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oblique travelling wave solution at finite amplitude. We now perturb the wavenumbers as 




a = a 0 + 6Ai(X, Z,T ), 

(109) 



II 

OCX 

+ 6A 2 (X, Z,T), 

(110) 

and hence transform 

the partial 

derivatives in the following manner 



d 

d 

J . d 8 \ 



dx 

a °d9 + d \ Al de + dx) ’ 

(111) 


d 

„ d 

./ 8 8 \ 



dz 

fio d6 + b { A2 d6 + dz J’ 

(112) 


d 

„ d 

c ( _ d 8 \ 



dt 

1 ^ 

0 

sO 

1 

+ A Ol d0 + dr) + ‘ 

(113) 

We also expand the velocity field 

and pressure function as 


34= (1 -y 2 ,0,0)' 

r r% 

r + Uo(0, 

,y) + Su 1 (X,Z,T,e,y) + ..., 

(114) 

p = 

[-Te X + <- 

-1 (X,Z) 

1 + \po(Q,y) + qo(X, Z, 7 1 )] + . . . , 

(115) 


where q x are the additional pressure terms required to satisfy a constant mass flux through the 
channel. 

Now as before we can expand the Navier-Stokes equations in terms of the slow scale pa- 
rameter, giving a leading order system 


1 f) 

— flo^O 0 + <^0^0^06 + VoUoy + PoWqUqq — V 3 U 0 + OLoPoB = ^ ~ , 

Re aX 

— floVoe + cxqUqVqb + vovoy + fioWo Voe — — V 3 U 0 + p 0y = 0 , 

He 

1 f) 

-tt 0 wo$ + aoUoWoe + VoWq v + fioWoWo# — — V^o + ftoPoe — , 

Re dZ 

OLqUqq + Voy + floWQQ = 0 , 


(116) 

(117) 

(118) 
(119) 



(«2 + /?o) 


d 2 


+ 


d 


do 2 dy 


( 120 ) 


where uq U -\-uq. At this order we let < 7—1 (AT, Z ) — K\X k 2 Z and choose /{j, k 2 to satisfy zero 
mass flux through the channel for the disturbance. Continuing further yields an 0(6) system 
of the form 


L u {u u v u w-i,pi} = Ai 
+ A 2 

L v {u 1 ,v 1 ,w 1 ,p 1 } = Ai 


—UoUoe — p 0 g -)- — — aoUogg 

He 


—wqUqs -f -=—f3oUoe& 
He 

2 

—uqVqq + — — oiQVoee 
He 


+ fliUoe — 


dqo 
dX ’ 


( 121 ) 
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r 2 

+ A 2 —wqvm + -^-fiovose + $1\Voq , 
L Re 

L w {u u v u wi,pi} = Ai 


+ A2 

# 0^10 + v ly + Pome = ^1 

where L u {ui, Vi,tui,pi} is given by 


—u o wo0 + —OL 0 wooe 
Re 


—wowoe — Poe + ~^r fio w oee 
Re 


+ ttiWoe — 


dqo 

dz 


-me 


+ A ^ 


-w 0 $ 


L u {ui,V\,W\ , Pi } = — fluio + <*0^0^10 + mU\Uod + VoUjy + V^Uoy 

+ PoWoUie + fiommo — — V3U1 + otoPie ? 

tie 

with analogous forms for L v and L w . 

After some work we can see that the solution is of the form 


(«i, t>i , m , Pi ) T = Aj ( X , Z, T) + ^ 12 A 2 (X, Z, T) 


ggo . 


+*,^(*,z,r) + * 1< ||(x,z,T), 


and the solvability condition determines fl\ as 


n, = —MX, z, T) + ^ a 2 (x, z, r) ; 

aao tfpo 


where we have used 


_ e&o 9 ki , ok 2 

Mi 


- dK *i 

- 7 T— ^ 


3 a 0 da 0 ~ 13 5 a 0 - 14 

_ ^<£0 < 9*1 d « 2 

' 10 Ck /O 0,0 


- 12 c«„ a *- 13 

»&l 

{ao,/?o^o} fixed, 


^,= 


^ 3 &C, 

~ 14 5 /C 2 


{ao,/?o,fio} fixed, 


( 122 ) 

(123) 

(124) 


(125) 


(126) 

(127) 

(128) 

(129) 

(130) 

(131) 

(132) 


with <f> 0 replaced by {tio, no, wo, Po}- 

In the usual two dimensional iteration scheme (Re fixed) we specify {ao, Amplitude} and 
iterate on {flo, /ci }■ Now for given wavenumber and amplitude, we can think of flo as a function 
of /cj with the required frequency given by f!o(/c{), where acj = /Cj is determined by the zero 
mass flux condition for the disturbance. Obviously we could equally fix {a 0 ,fIo} and iterate 
on {Amplitude, /ci} to define the amplitude as a function of the pressure constant. We can 
apply this same argument to the three dimensional iteration scheme and so define the solutions 

{£i3>£u)- 


24 


If we continue the expansion scheme to next order we obtain a similar system of equations 
with inhomogeneous terms of the form {A?, A 2 , A 1X , A 2X , A l9o *, . . we state the results for 
the fully three-dimensional problem in Appendix A but to simplify this problem (and outline 
the basic method) we shall consider the weakly three dimensional limit of 0 O -* 0. In this 
limit the leading order system will be the same two dimensional problem computed in previous 
sections but the solutions to the 0 (£) equations are 


where 


(u 1 ,v 1 ,u;,,p 1 ) T = (u„,n n ,0,p 11 ) T A 1 + (0, 0, w 12 , 0) r A 

' dq o 

dX ' v -’ dz 


+ (“ 13 , “ 13 , o, Pi3) T 7^7 + (0, 0, w 14 , 0) T |^ , 


(133) 

(134) 


(“n,“n,Pn) T = 

(“ 13 , “ 13 ,P 13 ) T = 


W 12 = 


“>14 = 


^-(u 0 ,n 0 ,po) r 
dw 0 8 k 2 

Wo~Wo WU ’ 

dw 0 | 


dtti 

da 0 


(“i 3 ,“i 3 ,Pi 3 ) r , 


{a 0 ,Po.si 0 } fixed, 


9k 2 


{a 0 ,Po,fio} fixed. 


The solvability condition at this order reduces to 


(135) 

(136) 

(137) 

(138) 


fti = |^A r(X,Z,T) 

OOCq 


(139) 


since for 0 O ~ 0 ( 7 ), 7 < 1 , 


0 ^ 

^-{u 0 , “0,Po, «l}, {Pl3, “ 13 , “13} ~ 0(1) , ^-{u>o, K 2 }, “>13 ~ 0 ( 7 ) , 

Q 

{“0,“0,P0, «l}, {Pl 4 , “14, “14} ~ 0( 7) , -^-{“> 0 , K 2}, “>14 ~ 0(1) . 


d0 o 


90 0 


(140) 


The group velocity obtained in this limit is equivalent to the two-dimensional case, as ex- 
pected. Continuing the expansion scheme further, to 0(8 2 ), yields the following after some 
simplification, 


L u {u 2 ,v 2 ,p 2 } = 


9Q 0 


L da Q Une ~ ™ oUne ~ — unuoe 

— v n u iiy — Pi 19 + — (2a o uii00 + u 0 ^)J Aj 


+ 

+ 


— W12U0Q + 


Re 



A 


2 

2 


-cxoU^u X3e - Vi 3 u X3y 


(9qo\ 2 

- 

Hax) + 

-^13^ 


^2 


9q 0 
9Z 
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-u 13 0 — u o ui 3 0 — a 0 unUi 3 e — a 0 «i 3 “n« — “13^00 


VllUl3y — Vl3U U y — Pi 3$ + — 


A,— 

l dX 


\dtlo . 1 , 9 

+ [ da 0 Ul3 U ° Ul3 Pl3+ R e ( 2a ° Ul30 ) IdX 2 

dSl 0 . 1 J 0Ax 

+ ^-«n — uo«n “ P11 + — [ZocoUne + u o 0 jj - qj £ 

, f . 1^2 dq 1 

+ [ /?e U H dZ + 2U(W ax’ 


(141) 


with an analogous form for the u-momentum equation and a decoupled iy-momentura problem. 
Continuity of mass at this order gives the further equation 

dv-2 dv 2 dAi d 2 qo dA 2 d 2 qo 

a °- Se + air = “ ““ Tx~ “ ,3 aF " Wu ~dz ~ W "W 


- U 11- 


d 2 q 0 dA 2 d 2 qo 

u ' 3 dT 2 ■ Wl2 az ■ WlA ~dz 2 


’u\ieA\ — U]3^A] ■ 


Wl 2 0A\ — Wi A gA 2 —^r . 


(142) 


We have used L u (in 141) to indicate the operator L u , defined by (125), with (3 0 = 0, and 
the boundary conditions for the systems discussed above are no-slip and impermeability at the 
fixed planes y = ±1. In this three dimensional formulation the additional pressure term, qo, 
will contribute explicitly to the solvability condition for S7 2 ; whereas in the previous theory this 
term was formally eliminated using the stream function and vorticity equation approach. This 
is analogous to the work of Davey, Hocking and Stewartson [4], concerning the weakly nonlinear 
evolution of three dimensional disturbances in plane Poiseuille flow, where they found that a 
secular pressure term contributed to the final form of the evolution equation. 

If we consider the continuity equation for those terms which have zero mean with respect 
to the phase variable 0 , we obtain 

dv 2 _ dAi _ d 2 q 0 _ dA 2 _ d 2 q 0 . ^ 

~d^~~ Uu dX Ul3 dX 2 Wu dZ ™ XA dZ 2 ' ^ ’ 


_ d 2 q 0 
U ' 3 dX 2 


(143) 


where the bar notation indicates the mean part of the relevant expression. Integration of this 
equation shows that, for impermeability at both boundaries, we must satisfy the linear partial 
differential equation 

d 2 qo T d 2 q 0f d Aid hi dA 2 dn 2 , 


^ T , v W T _ v * 

dx 2 1 dZ 2 2 “ dX da c 


where 


dA 2 dn 2 
+ dZ d/3 o' 

(144) 

/*+! 

/ w 14 d y . 

Jy =- 1 

(145) 


Now returning to the solvability condition determining H 2 , we can see that further simpli- 
fications are possible by introducing some particular solutions; for example, 

(>46) 

l UKi OK\ UK>\ J 


(146) 
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(147) 

(148) 


L-ll, 


dv\3 dp\3 \ _ „ 
8ki ’ dni ’ 3ki J 2 ’ 
duu dvi2 dpu\ _ c, 

3k 2 ’ 3k 2 ’ d/c 2 f 3 ’ 


where Si,S2, 53 are the [. . .] coefficients of the (<7o;r) 2 an d A 2 %z terms respectively in 

(141). These particular solutions can be applied to the system as a whole, reducing the form 
of the frequency correction term to 


Cl 2 — f^2i Aj T f^A^ + 4*i 


8 Aj 

3X 


+ <F 


3 A 2 
: 3Z 


+ 4*3 


d 2 <? o 

3X 2 ’ 


(149) 


where {fl 2 !, 4>i} are constants determined from the appropriate integrals contained in the solv- 
ability condition. Furthermore a perturbation of {a 0 , A)} in the leading order problem shows 
that 

n _ 1 d 2 Q 0 ^ 1 3 2 Q o 

21 "2 3a 2 ’ ° 22 “ 2 3& 2 ‘ 

The consistency conditions now determine the evolution of the wavenumber perturbations, and 
applying the above results we obtain 



3A^ dSl 0 3Ai 2 

S ~dT +S d^-dX +S 


It"' 

3 2 Ai 


3A 1 3 2 n o 3A 2 

+ A 2 - 


3X ' 381 
3 2 A 2 


3X 


^ 3 2 A ! ^ <9 2 A 2 ^ a 3 9 ol _ _ 

1 3X 2 2 3X3Z^ 3 3X 3 )^ { ’ ’ 


3Ai _ 3A 2 
3Z ~ ~3X ’ 


(151) 

(152) 


Now for a leading order balance we again introduce the slower timescale r (defined by 40) and 
make the Galilean transformation (42). Now by defining 0 as 


<90 

di 


= A a 


<90 

3Z 


= A 


2 , 


we obtain 


d 2 fl 0 - - <9 2 f> 0 - - 

+ + ~Qp2 = -4>10^ — 4> 2 ®(ZZ - 4 > 3 <? 0 £(( , 


(153) 


(154) 


and 


+ Siq 0 zz = — A 2 0^ — A 3 0^z • 


(155) 


As mentioned above the terms {4>i, 4> 2 , 4> 3 } can be calculated numerically (although it is not 
a trivial computation) from the integrals derived in the solvability condition, and {A!,A 2 ,A 3 } 
denote the coefficients {/ 2 / /i, <9/Ci/da 0 , (/ 2 / h) 3 k 2 /38o} in (144). We observe that if 8z = 0 
then qozzz reduces to a constant multiple of Ai^, the O(^) problem is simplified, and we return 
to Burgers equation once more. 
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There is little we can say about the solution of this equation without a detailed numerical 
investigation and computation of the coefficients from their integral representations. The form 
of the system (154)— (155) can be simplified by rescaling appropriately but only to the form 

A X t + A X A XX ± A(A( X iA xxx + c-iA x ^ + B XXX ’ (156) 

B xx A C-2 -^CC — ^2-4xx T £zA^ , 0^^) 

where Cj are constants, which can be written in terms of those appearing previously. We 
note that the first equation, (156), can be integrated with respect to x if we redefine the 
dependent variable to include an arbitrary function of r. We also note (as in Davey, Hocking 
and Stewartson [4]) that a return to Burgers equation (after a further rescaling) is achieved 
for skewed two-dimensional perturbations, i.e. for wavenumber perturbations which may be 
written solely in terms of a skewed variable 

E = /x + m(, (158) 

for some real values / and m. Obviously the sign of the diffusive term will determine the type 
of solution obtained, but this must again be calculated from the coefficients introduced above 
and any chosen values of {/,m}. 


7 Conclusions 

We have applied a phase equation technique to develop a perturbation theory for slowly varying 
finite amplitude wavetrains, at finite Reynolds number, in plane Poiseuille flow. When the 
method is applied to uniform wavetrains we have shown that a small wavenumber perturbation 
evolves according to Burgers equation. We can obtain an exact solution to Burgers equation 
via the Cole-Hopf transformation and the solution is characterized by the sign of the diffusive 
coefficient. For a typical uniform wavetrain with 0(1) disturbance energy we have discussed 
how to compute the coefficients appearing in the evolution equation. 

A number of numerical results, for solutions corresponding to cross sections of the neutral 
surface, have been presented showing that both diffusively stable and unstable cases are pos- 
sible in general. In the unstable case ($ positive) we know that the wavenumber perturbation 
develops a singularity at finite time in the slow scale. For the diffusively stable case ($ nega- 
tive) it is possible for weak shock structures to appear for non-localized or non-periodic initial 
conditions. We make no claims that this instability of the finite amplitude wavetrain is in any 
way the “most unstable'”, but since an initial disturbance can always be found that does not 
decay, it suggests that the use of uniform wavetrains in theory and computation needs some 
justification for large scale problems. 

As the shock structures in the wavenumber perturbation are approached we expect higher 
order spatial derivatives to become of equal importance in the evolution equation. We see (as 
noted by Bernoff [1]), by an induction argument, that higher order frequency corrections are of 
the form 

Q — P 

— 1 n 


da d 2 a 
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dX' dX 2 '"’' dX n 
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a polynomial in the slow scale derivatives of the wavenumber. Thus as the length scales shorten 
all the higher order terms will become of comparable magnitude simultaneously. When all the 
previously neglected terms return to the evolution equation we must return to the full unsteady 
2-D Navier-Stokes equations to determine the development. 

We have shown in §3 that the phase equation method is distinct from a multi-scale amplitude 
perturbation approach. A perturbation of the phase in the form 

9 = 9 Q (x,t) + ^® x (8x,8t) + ... , (160) 

may be expected to be equivalent to an amplitude perturbation for e sufficiently small. For 
uniform wavetrains however we showed, in §2, that for a leading order balance we required that 
the wavenumber perturbation be of comparable magnitude to the slow scale parameter, thus 
e ~ 8 and the multi-scale approach must be distinct. The evolution equation obtained from the 
amplitude perturbation method is simply the heat transfer equation and no nonlinearity can 
be introduced. For the diffusively unstable case we must rely on linear higher order derivatives 
to become of comparable magnitude and alter the form of the evolution equation. 

We have also shown that the evolution equation must breakdown to a different form as 
the leading order problem approaches the linear limit or the region where the upper and lower 
branches join at finite amplitude. This singularity is due to the behaviour of the neutral 
surface in these areas, and in the weakly nonlinear case we observed that a regularization of the 
evolution equation is nontrivial. Indeed this has to be the case if we are to be able to reconcile 
the phase equation approach with the well known results of weakly nonlinear theory, namely 
the Stuart- Watson equation. 

To investigate the effects of three dimensionality we also developed a form for the evolution 
equation by extending the definitions of the phase method applied in previous sections. To 
simplify the basic analysis a weak spanwise dependence is allowed for in the uniform wavetrain 
problem; the analogous results for a fully three-dimensional situation are stated in the Ap- 
pendix. In this case we show that the evolution equation is more complex, and is coupled with 
a further linear partial differential equation. This coupled equation determines an additional 
pressure function (required in order satisfy an impermeability at the parallel boundaries) that 
in the two-dimensional formulation is formally eliminated using a stream function approach. 
The comments made above concerning the breakdown of the phase equation method, at the 
E = 0 plane of the neutral surface, will similarly apply here in the three-dimensional approach. 
We should expect to be able to, likewise, match the phase method, in the appropriate areas of 
the neutral surface, to the results of Davey, Hocking and Stewartson [4] concerning the weakly 
nonlinear development of three-dimensional disturbances. The multi-scale approach, of §3, can 
be similarly extended to three dimensions, giving the heat equation with diffusion on both slow 
scales although this is again linearly coupled with the same partial differential equation for the 
pressure term. 

Thus generally the development of, and role played by, these uniform wavetrains in PPF is 
complex to determine. We have shown susceptibility to weak shocks and finite time singularities 
for both upper and lower branch solutions, together with a complex regularization problem in 
the E — > 0 limit. As well as these slow scale effects we also have the superharmonic instability 
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results of Pugh k SafFman [14], which occur on 0(1) scales, the apparent slow decay of three 
dimensional perturbations to the fully developed flow as discussed by Orzsag k Patera [13] and 
the bypass mechanisms of Gustavsson et al , to mention just some of the instability mechanisms 
present. Obviously it is difficult to say (in any given problem) how such mechanisms will 
interact but to assume that any one will dominate over others needs to be carefully justified. 

There are a number of questions that remain unanswered in this discussion; i.e. matching 
the phase equation theory to a relevant approach in the weakly nonlinear case, and the effect 
of the weak three dimensionality through computation of the evolution equations (154-155). If 
in future we wish to produce a full numerical procedure, to give refined values for the neutral 
surface, or for computation of the evolution equation with weak spanwise dependence, we should 
perhaps consider replacing the numerical method described previously with a similar collocation 
method. Such techniques were used successfully at a later stage by Herbert [9] and have the 
advantage of reducing computation time by allowing for a simpler evaluation of the nonlinear 
terms in Newtons method. 

The authors would like to acknowledge the financial support of EPSRC for R.E.H. 


A The Fully 3— D Problem 


In this appendix we present some of the more lengthy details involved with the application 
of phase equation methods to a fully three-dimensional problem; by which we mean (in this 
application) the modulation over long spanwise and streamwise scales of a fully oblique uniform 
wavetrain in PPF. We prefer to present these details as an appendix in order to keep the basic 
method (as presented in §6 for the case of a weak spanwise dependence) as simple as possible. 

Substituting the form of solution for the 0(6) system, defined by (129-132), into the inho- 
mogeneous terms at 0(6 2 ) yields 
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as the vector form of the momentum equations, where L = (L U ,L V ,L W ) T , 
continuity of mass condition 
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The coefficients, denoted by R, above, are straightforward to determine; for example 
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This system requires a solvability condition to be satisfied and therefore determines fl 2 for a 
given leading order solution. Before computing the integrals needed for the evaluation of fl 2 
we note that a number of the inhomogeneous terms given above can be removed from the 
solvability condition by introducing some particular solutions, 
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together with the previous solutions (129)— (132); obviously these solutions likewise apply to 
the appropriate inhomogeneous terms of the continuity equation. 

The explicit content of the solvability condition can be determined by computing the adjoint 
and relevant integrals, although we do not give such details here it is noted that the final form 
of the frequency correction must be 
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for coefficients fl 2l and 4> t , which in general must be determined numerically. Although we can 
not remove the terms A 2 , A]A 2 , A^ from the solvability condition we can effectively remove 
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them from the computation by noting that a further wavenumber perturbation of the systems 
satisfied by {wn, wn> an< ^ { u i 2 ; ^ 12 , ^ 12 ,^ 12 } yields 
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The evolution of the wavenumber perturbations is then governed by the following system 
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where we have introduced 
defined by 


the constants A,-, the slower timescale r 
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and the new coordinate ( 
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The equation (172) and constants A; are obtained from solving the continuity equation (162) for 
the mean flow term, v 2 , and applying the impermeability condition at the boundaries y = ±1. 


References 

[1] Andrew J Bernoff. Slowly Varying Fully Nonlinear Wavetrains in the Ginzburg-Landau 
Equation. Physica D 30 1988 

[2] K.M. Butler, & B.F. Farrell. Three Dimensional Optimal Perturbations in Viscous Shear 
Flow. Phys . Fluids A . 4 1992 

[3] M.C. Cross k A.C. Newell Convection Patterns in Large Aspect Ratio Systems. Physica D 
10 1984 

[4] A. Davey, L.M. Hocking k K.Stewartson. On the Nonlinear Evolution of Three-Dimensional 
Disturbances in Plane Poiseuille Flow. J. Fluid Mcch. 63 1974 

[5] L.H. Gustavsson. Growth of Three Dimensional Disturbances in Plane Poiseuille Flow. J . 
Fluid Mech. 224 1991 


32 



[6] L.H. Gustavsson. Resonant Growth of Three Dimensional Disturbances in Plane Poiseuille 
Flow. J. Fluid Mech. 112 1981 

[7] L.H. Gustavsson & L.S. Hultgren. A Resonance Mechanism in Plane Couette Flow. J Fluid 
Mech. 98 1980 

[8] Philip Hall. A Phase Equation Approach to Boundary Layer Instability Theory: Tollmien 
Schlichting Waves. J. Fluid Mech. to appear 

[9] Thorwald Herbert. Periodic Secondary Motions in a Plane Channel. Springer Lecture Notes 
in Physics 59 1976 

[10] L.N. Howard and N.Kopell. Slowly Varying Waves and Shock Structures in Reaction- 
Diffusion Equations. Stud, in App. Math. 56 1977 

[11] A. C. Newell, T. Passot. Towards a Universal Theory for Natural Patterns. Physica D 74 
1994 

[12] A.C. Newell, T. Passot, J. Lega. Order Parameter Equations for Patterns. Ann. Rev. Fluid 
Dynamics 25 1993 

[13] S.A. Orzsag , A.T. Patera. Three Dimensional Instability of Plane Channel Flows. J. Fluid 
Mech. 128 1983 

[14] J.D.Pugh and P.G.Saffman. Two-Dimensional Superharmonic Stability of Finite- 
Amplitude Waves in Plane Poiseuille Flow. J. Fluid Mech. 194 1988 

[15] L.N. Trefethen, A.E. Trefethen, Sz S.C. Reddy. Pseudospectra of the Linear Navier-Stokes 
Evolution Operator and Instability of Plane Poiseuille and Couette flows. Dept. Computer 
Science, Cornell University, preliminary report. 1992 

[16] G.B.Whitham. Linear and Nonlinear Waves. Wiley- Interscience, New York 1974 

[17] J.P. Zahn, J. Toomre, E.A. Spiegel, D.O. Gough. Nonlinear Cellular Motions in Plane 
Poiseuille Channel Flow. J. Fluid Mech. 64 1974 


33 



Form Approved 
OMB No 0704-0188 


REPORT DOCUMENTATION PAGE 


Public reporting burden for this collection of information is estimated to average 1 hour per response, mcludmgthe time for reviewing instructions, searching existing data sources. 
Mtherinjrand maintaining the data needed, and completing and reviewing the collection of information, Send comments regarding this burden estimate or any oth« r ViCUffiilon 
rnll^rtmn of information incl udine suegestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports 1215 Jefferson 
Davis Highway, Suite 1204, Arlington?^* 22202-4302, and to the Office of Management and Budget, Pap erwork Reduction Project (0704-0188), Washington, DC 20503. 


1. AGENCY USE ONVf(Leave blank) I 2, REPORT DATE 


3. REPORT TYPE AND DATES COVERED 

Contractor Report 


«. TITLE AND SUBTITLE 

THE EVOLUTION OF FINITE AMPLITUDE WAVETRAINS IN 
PLANE CHANNEL FLOW 

5. FUNDING NUMBERS 

C NAS1-19480 
WU 505-90-52-01 

6. AUTHOR(S) 

R. E. Hewitt 
P. Hall 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Institute for Computer Applications in Science and Engineering 
Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23681-0001 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

ICASE Report No. 96-10 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681-0001 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA CR-198281 
ICASE Report No. 96-10 

11. SUPPLEMENTARY NOTES 

Langley Technical Monitor: Dennis M. Bushnell 
Final Report 

Submitted to Phil. Transactions of the Royal Society 

12a. DISTRIBUTION/AVAILABILITY STATEMENT 

U nclassified-U nlimited 

12b. DISTRIBUTION CODE 

Subject Category 34 



13. ABSTRACT (Maximum 200 words) 

We consider a viscous incompressible fluid flow driven between two parallel plates by a constant pressure gradient. 
The flow is at a finite Reynolds number, with an 0(1) disturbance in the form of a traveling wave. A phase equation 
approach is used to discuss the evolution of slowly varying fully nonlinear two dimensional wavetrains. We consider 
uniform wavetrains in detail, showing that the development of a wavenumber perturbation is governed by Burgers 
equation in most cases. The wavenumber perturbation theory, constructed using the phase equation approach for a 
uniform wavetrain, is shown to be distinct from an amplitude perturbation expansion about the periodic flow. In 
fact we show that the amplitude equation contains only linear terms and is simply the heat equation. We review, 
briefly, the well known dynamics of Burgers equation, which imply that both shock structures and finite time 
singularities of the wavenumber perturbation can occur with respect to the slow scales. Numerical computations 
have been performed to identify areas of the (wavenumber, Reynolds number, energy) neutral surface for which each 
of these possibilities can occur. We note that the evolution equations will breakdown under certain circumstances, 
in particular for a weakly nonlinear secondary flow. Finally we extend the theory to three dimensions and discuss 
the limit of a weak spanwise dependence for uniform wavetrains, showing that two functions are required to describe 
the evolution. These unknowns are a phase and a pressure function which satisfy a pair of linearly coupled partial 
differential equations. The results obtained from applying the same analysis to the fully three dimensional problem 
are included as an appendix. 


14. SUBJECT TERMS 

wavetrains; instability; channel flows 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

NSN 7540-01-280-5500 


15. NUMBER OF PAGES 

35 


16. PRICE CODE 

A03 


14. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION 
OF THIS PAGE OF ABSTRACT OF ABSTRACT 

Unclassified 

Standard Form 298(Rev. 2-89) 

Prescribed by ANSI Std, Z39-18 
298-102 












